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In this methodology focused paper we scrutinize the application of the band-coherent Boltzmann 

equation approach to calculating the conductivity of chiral particles. As the ideal testing ground 

}_^ we use the two-band kinetic Hamiltonian with an TV-fold chiral twist that arise in a low-energy 

<**> | description of charge carriers in rhombohedrally stacked multilayer graphene. To understand the 

■^^ role of chirality in the conductivity of such particles we also consider the artificial model with the 

chiral winding number decoupled from the power of the dispersion. We first utilize the approximate 
but analytically solvable band-coherent Boltzmann approach including the ill-understood principal 
value terms that are a byproduct of several quantum-many body theory derivations of Boltzmann 
collision integrals. Further on, we employ the finite-size Kubo formula with the exact diagonalization 

,—1 of the total Hamiltonian perturbed by disorder. Finally, we compare several choices of Ansatz in the 

derivation of the Boltzmann equation according to the qualitative agreement between the Boltzmann 
and Kubo conductivities. We find that the best agreement can be reached in the approach where 
the principle value terms in the collision integral are absent. 



I. INTRODUCTION 



The isolation of the one-atomic carbon layers [T] has boosted an extremely intensive research in the field of quantum 
transport associated with Dirac electrons. [5] One of the primary reasons of such an overwhelmed interest is the bizarre 
chiral nature of the Dirac carriers which makes the conduction and valence states entangled with each other. This 
C^ coupling leads to many nontrivial effects, [2 including the famous Klein tunneling 3J. 

The chiral carriers are usually characterized by means of the pseudospin which is an additional quantum number 

i ^-*, formally similar to the real spin in spin-orbit coupled systems. In a single layer graphene, the pseudospin keeps its 
radial orientation along the Fermi circle resulting in the winding number N = 1. The winding number is directly 

t-H related to the Berry phase irN responsible for the anomalous quantum Hall effect observed at the very beginning of the 

graphene rush. [H [5] There are quite a few speculations on the important role of the pseudospin-coherent contribution 
in the electrical dc conductivity of Dirac fermions at low carrier concentrations. [rjHTU] There is however no agreement 

i^J whether this contribution really matters and can be extracted from the total conductivity. 

-y~. There are two aspects of this problem worthy of discussion. The first one is experimental: It is very difficult experi- 

mentally to detect the pseudospin-coherent contribution because the subtle quantum mechanical effects responsible for 
this contribution are easily smeared out by temperature, ripples, and charged inhomogeneities inherited in graphene. 
In spite of the impressive progressjll) made with suspended graphene, this contribution has not been extracted yet. 
The second problem is purely theoretical: It is still not clear how large the pseudospin contribution should be since 
• • the numerical methods can give only an estimation [T^] while the approximated Boltzmann-like approaches [SHIP] are 
unreliable at low carrier concentrations and may lead to different answers depending on the approximations involved. 
The main goal of the present work is to compare the exact numerical (finite-size Kubo) and approximated analytical 
(Boltzmann-like) pseudospin-coherent conductivities and in that way to figure out the correct approximation needed 
to derive the most reasonable analytical expression. In particular, we focus on the dependence of pseudospin-coherent 
conductivity on the winding number N. Surprisingly, the semiclassical pseudospin-coherent conductivity acquires 
qualitatively different behavior on N depending on the approximation done. The subsequent comparison with the 
numerical calculation helps us to rule out the "wrong" approximations and keep the "correct" ones. 

Note, that the chiral particles with an arbitrary winding number can be realized in the chirally stacked multilayer 
graphene. The recent investigations of the band structure in a few layer graphene [13H15] indicate its strong sensitivity 
to the stacking pattern. Graphene layers placed together in the most natural way do not lie exactly one on top of 
the other but are shifted in such a way that only half of the carbon atoms have a neighbor in another layer and the 
other half are projected right into the middle of the graphene's honeycomb cell. If the third layer aligns with the 
first (and the N + 2 layer with the TV-th) then we arrive at the most abundantly occurring arrangement of graphene 



layers known as Bernal (or ABAB...) stacking. An alternative stacking is obtained when the third layer aligns neither 
with the first nor with the second but is shifted with respect to both, and it is the fourth layer that again aligns with 
the first. This arrangement is known as rhombohedral (or ABCABC.) stacking. The simplest model used in the 
study of the latter system [TB] is the model that we shall use as our testing ground in the investigation on theoretical 
methods to calculate the conductivity. 

The beauty of the effective low-energy Hamiltonian describing the carriers in the chirally stacked multilayer graphcne 
has attracted many theoreticians working in the field of quantum transport. Among the most recent works are 
the theory of chirality-dependent phonon scattering [17) . the semiclassical Boltzmann transport theory |18j . the 
investigations of Landau levels [15] ■ Berry phase [20], effects of screening J3TJ[22], optical conductivity [2"3~] . and THz 
absorption [24]. We would like to warn the reader, that although our ultimate interest is to understand conductivity 
in graphene type of systems, the aim of this paper is purely methodological, and we undertake no discussion about 
the connection to realistic systems and to experimental findings. For example, we will for simplicity and analytical 
solvability assume a fully homogenous background, whereas realistic graphene might have to be modeled with a 
variable potential to account for ripples and other sources of inhomogeneity. 

The outline of the paper is the following. In section ITT] we discuss some preliminaries that will be used in later 



sections. Section IV offers in more detail some of the background and motivation of the present paper. In section III 
we present an analytical heuristic derivation of the conductivity by evaluating the Kubo formula using the chiral 
plane-wave states of free electrons of multilayer graphene. The derivation, interesting in its own right, will serve 
among other as a easy introduction to band-coherent effects and to the distinction between intra-band and inter- 
band contributions to the conductivity. In section [V] we will discuss the conductivity derived with a band-coherent 
Boltzmann equation approach. In section |VI| we present the results of the numerical studies of the Kubo formula. In 
section [VlTJ we wrap up with a discussion. 

II. PRELIMINARIES 

Many properties of bilayer graphene can be well understood in terms of two independent, time- reversal related copies 
of the effective two-band Hamiltonian [2"5] 

rr _ l ( (p x -iPy) 2 \ m 

This hamiltonian exhibits a coupling between momentum and a spin type of variable, a pseudo-spin that in this case 
derives from the non-equivalent lattice sites, one from each layer. For ./V-layer rhombohcdrally stacked graphene one 
can with similar approximations derive an effective low-energy hamiltonian that also turns also out to be a two band 
model, this time given by [T|)] 

H °- 1 [( Px + lP y) N J' (2) 

with 7 (like m* above) being determined by material parameters of the system. As a realistic models for rhombohedral 
multilayer graphene the hamiltonian (pi is expected to become increasingly poor with an increasing number of layers. 
We wish to exploit the hamiltonian (pi for theoretical and methodological studies because of its combination of 
simplicity and yet very non-trivial chiral properties. For our studies it will turn out to be interesting to generalize 
the hamiltonian (pi to the more general although artificial hamiltonian 

H = b<r = b(k)b(k) ■ a = j(hk) N * ( e X e e "^° 9 ) , (3) 

where fc = |fe| is the absolute value of the wave number k = kk = p/H, 8 is the angle of the direction k — k/k — 
(cos 8, sin 8) of the wave number and 6(fe) is the (pseudo)spin-orbit field. In the hamiltonian (pi) the power N& of k in 
the dispersion e& = b = r y(hk) Nd has been decoupled from the winding number N c of the direction b of the pseudospin- 
orbit field as a function of k, which determines the chiral properties of the hamiltonian. In the special case N$ = N c 
this hamiltonian reduces to the hamiltonian (pi . The benefit of introducing such an artificial Hamiltonian is that one 
can take advantage of some special features of a quadratic dispersion, to be explained later, but still keep the chiral 
properties general. 

For both the numerical and analytical treatments of the conductivity we will need the velocity matrix, defined as 
Vj = 'W 2 -- For the multilayer Hamiltonian (pi one has in the basis of the sublattice orbitals (that is, in the basis in 



which equation (pi) is written) the x-component 

For us it will sometimes be convenient to express matrix quantities in the eigenbasis of the Hamiltonian, which we 
also call the chiral basis. For the Hamiltonian fcfy the eigenstates are given by 



'ik-x 



where s = ± is the pseudospin index and where 9 = aictani{ky / k x ) is the angle of the momentum vector in the 
momentum plane. The eigenenergies are e| = sb = sj(hk) Nd . The velocity operator in the chiral basis reads 

k x d k b iO x ^Y h _ ( cos9 _^sin( X " 



Vx = k x h- 1 d k Ho + e x (hk)- 1 deH = h- 1 [ ™; w 2 fc = «* .^ . "- (6) 
v ' s 1 ' / V -iQ^ ~kdhb I \ lysine* -cos9 / 
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where 9 = z x k = (— sin#,cos0) and v k = d k b = N < x , y(Hk) Nd ~ l . The velocity operator can formally be decomposed 
into two parts. The intraband part v^. ntra accounts for the ordinary velocity components for each band as where 
the bands not talking to each other. In the chiral basis it is diagonal. On the other hand, the interband part 
v™ ter = 0(hk)^ 1 deH couples the two bands and is an off-diagonal matrix in the chiral basis. The interpretation 
of these matrix components is less intuitive. They reflect the Zitterbewegung property of Dirac electrons. Both 
intraband and interband components of the velocity can contribute to the current. Within the Kubo formula picture 
we will see this in next section. In the Boltzmann picture, we note that for a matrix-valued distribution function / 
the current in the ^-direction is given by 

3x = J ^Tr (y w f) - 2 J v k (k x f B + 6 x h^\ (7) 

in terms of the components of the matrix- valued distribution function 

/ = /ol + / • <r = /ol + (/&& + he + /,£)• <r = ( J* + ,-| ff~_ if f * ) (8) 

with the vector part or spin part of / decomposed along the ^-dependent basis vectors {b, c, z} with c = b x z. The 
densities n = /o ± ft (in phase space) in the upper and lower band, respectively, we find in the chiral basis on the 
diagonal. Note however, that for the type of Hamiltonians we consider, only the spin component / but not the total 
charge density 2/o = n + + vT contributes to the current. 

The conductivities a — j x /E x (the driving electric field is assumed to be in the ^-direction) will in the present 
paper be given as functions of the dimensionless variable 

£kp = T tI vpkp = Nd — - — . (9) 

n 

The momentum relaxation time r tr = (T(k-p))^ 1 is derived from integral 

/d 2 k' N A(9 

j^yS(e k e k ,)W kk , cos 2 -^(1 - cos A9) (10) 

where W kk > = ^ni\U kk 'L 2 \ 2 gives the collision probability rate according to Fermi's golden rule for spinless particles 
scattered by an impurity potential U and the unusual trigonometric factor cos 2 —\ — comes from the spin overlap 
between the chiral eigenstates at a momentum angle difference A9 = 9 — 9' . In particular, the Drude conductivity is 
given by 

e 2 ik F 

ffD = X'^' (11) 



independently of the dispersion and winding number. 



The Kubo formula for the conductivity reads 



.fie 2 -r^ /fd (en) - /fd (e n ') (n\v x \n')(n'\v x \n) 
o-Kubo = -»tt > _, ;— , (12) 
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where \n) are the eigenstates of the total Hamiltonian including impurities but excluding the driving electric potential 
and /fd (e) is the Fermi-Dirac distribution. The parameter 77 describes the broadening of the eigenstates. For a pure 
sample it would be infinitessimal. In our case it is finite, as we describe a system with impurities. We will find it 
useful to study the quantities 

<7 lntra := CT Ku boK -> vj! ltra ) 

<x intor := CT Ku boK -> vf-) . (13) 

Numerical checks will show that there is at the most a negligible contribution from the remaining possible terms 
containing matrix elements (n|v! c ntra |ri')(n'|vj. ntcr |n) or (n|v™ ter |n')(n / |v™ tra |n) . 

It is rather clear that the Drude conductivity ctd as a pure intraband contribution resides completely within cr lntra . 
However, the above definitions of cr lntra and cr mtGr are in general not a clean decomposition of intraband and interband 
contributions to the conductivity, wherefore a lntra might on top of ctd contain some interband contributions. For 
the problem we study in the present paper the Boltzmann approach gives support for this decomposition faithfully 
separating pure intraband contributions form interband contributions, but in other cases, for example for the much 
more complex case of monolayer graphene (N =1), the Boltzmann approach (see ref. |H|) comes with the nonintuitive 
message that there can be important spin-coherence originated contributions also in the part /t of the distribution 
function that contributes to the current through v; ntra , although ft = n + — nT at a first glance seems to be a pure 
intraband contribution as it describes the polarization along the quantization axis of the Hamiltonian. 

III. HEURISTIC EVALUATION OF THE KUBO FORMULA USING THE PLANE WAVE BASIS 



In the analytical derivation based on the Kubo formula ( 12 ) presented in this section we will assume that a lntI>i 
contains all intraband contributions, in this case only the Drude conductivity ctd, whereas a mtcr contains all interband 
contributions. (As mentioned, support for this assumption will come from the Boltzmann approach.) We therefore 
choose r\ such that we obtain cr lntra = o- D and then derive <7 lntor with the given value of 77, hoping that we obtain a 
sensible result. This is one of the assumptions make the derivation a heuristic one. 

What furthermore allows an analytical treatment but also contributes to the heuristic nature of the derivation is 
that we assume that we — rather than using the eigenstates of the full problem including impurities — can use the the 
eigenstates W\ of the free Hamiltonian, that is the plane wave states, which makes an analytic evaluation simple. At 
high densities, that is, at high Fermi momenta, the use of plane wave states should be a good approximation as the 
weak impurities deform the eigenstates only slightly from the plane wave states. However, at small Fermi momenta, 
close to the Dirac point, one would expect the eigenstates to be substantially deformed by the impurities and we 
know of no solid argument for why the use of plane wave states would give sensible results also close to the Dirac 
point. Nonetheless, the outcome appears to be sensible. In any case, the derivation gives a nice introduction into the 
non-intuitive interband contributions to the conductivity. 

The calculation is done in the chiral basis. For the evaluation we use that 

(s, s') (ks\v x \k's')(k's'\v x \ks) 



(+1,+1) v i^kk'^k'k( cos ^) 2 intraband 

(-1,-1) Wfc<W<W-cos6>) 2 intraband (14) 

(+1,-1) v\ S hk ,5 k > k (-i^sin9)(i^shi0) interband 

(-!,+!) vlS kk >5 k , k {i^sm9)(-i^sme) interband. 



Thus, we can compactly summarize this as 

l — ss' 



fe 2 V?;2 /FD(eiQ-/FD(ei-') UJ (1 WCOB26Q/2 

L z *—i el — el el — e- - 1 - '" 



"/ 



For the intraband part s' — s we find the zero e| — e s k in the denominator. However, also the numerator contains a 
zero in /fd (ej) — /fd^ )• Since the arguments of the two Fermi-Dirac functions are 'close' to each other, we can 



Taylor expand the difference as 



/FD(eD-/FD(eI) = + (4-4) 



d/FD (eg) 
del 



(16) 



We see that the numerator contains the same zero as the denominator and we can cancel the two occurrences of 



el to obtain the intralayer contribution 
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/ ttt^ cos ^ £ f - e) = / — cos 2 / deD(e> 2 5(e F - e) = -^-£>f« f = -^- - *fvf , (17) 

J {2-k) 1 r\ J 2tt J 2rj 2ft i] 



where we used D{e) = ^—f— and where we have assumed zero temperature. We set r\ = h/r tI in order for a" 
reproduce the Drude conductivity od- 

We now turn to the interband contribution. We find 
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(18) 



With the previously assumed rj — h/r tr we obtain 



Winter 



Ah N d 



( 2e F r tl 
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V h 



The total conductivity is therefore 



The rewriting 
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(19) 



(20) 



(21) 



with the arcus tangens function taylor expanded in powers (Ik-p)^ 1 shows how an infinite series of terms of pseudospin 
originated quantum corrections diverging at ik-p = can sum up to something finite. In the Boltzmann regime £kp 3> 1 
the interband part vanishes and we are left with the Drude contribution. In the opposite limit £kp = the Drude 
conductivity vanishes and the interband contribution saturates at the finite value 



= 0) = 



h 8 N d 



(22) 



Note that the increasing the chiral winding number N c pushes up the conductivity at the Dirac point, whereas 

increasing the power N d of the dispersion pushes it down. For the multilayer Hamiltonian (pi (7V C = N d = N) the 

total conductivity has a minimum at the Dirac point £kp = 0, and the value of this minimum is proportional to the 

2 

number of layers N. In particular, for the bilayer case N — 2 this agrees with the result <x m ; n = 4-f in ref . [261 derived 

analytically with a spin-coherent Boltzmann type approach. In the case N c > N d the qualitative behavior is different. 

A minimum occurs at some non-zero £fc F whereas close enough to Ik-p = the conductivity shoots up to saturate at 

a local maximum, with the value given by (22). Examples are shown in Figs. [I] and p] 
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FIG. 1: The Boltzmann conductivity for different Ad and N c . The leftmost panel shows the case of Ad = 2 (quadratic 
dispersion), where all three Boltzmann approaches agree. The other three panels show the Boltzmann conductivity for Ad = 3 
with different choices of Ansatz, that is, with different outcomes for the principal value terms in the collision integral. The gray 
line corresponds to the Drude prediction. These results suggest that it is the increase chiral winding number A c that pushes up 
the conductivity close to the Dirac point as on increases N, whereas the increase in density of states at the Dirac point that 
comes with increasing Ad instead decreases the conductivity. Note that this observation applies for all considered choices of 
Ansatz . 
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FIG. 2: The analytical results for the total conductivity for JV(= Ad = N c ) equal to 2 (red), 3 (yellow), 4 (green), 6 (blue), and 
8 (magenta). The Drude prediction is shown by the gray line. The leftmost panel shows the results of the heuristic analytical 
treatment of the Kubo formula, and the subsequent three panels demonstrate the Boltzmann conductivities for a particular 
Ansatz. Note that the conductivity at the neutrality point £k-p = decreases with increasing A in the case of the Ansatz AA, 
but does not do so in the other considered cases. Note also that at £kp — the values with the Ansatz GKBA agree exactly 
with those with the heuristic treatment, but that for other £kp there is no agreement for A > 2. The further comparison with 
the numerically calculated conductivity will allow us to figure out, which Ansatz should be considered as the most proper one. 



IV. DETAILED INTRODUCTION TO THE BACKGROUND AND MOTIVATION OF THE PAPER 

This section complements the introduction with a more technical and detailed account of the background and mo- 
tivation of the present study as for the methodological issues concerning the use of band-coherent generalizations of 
Boltzmann equations to study graphene-related systems. 

The Boltzmann equation is a kinetic equation that is a widely used starting point for deriving the transport 
coefficients and other response functions, also for quantum systems, where it is used as a semiclassical approximation. 
The appeal of a Boltzmann type of approach is that it makes contact with a classically intuitive picture, in contrast 
with the kind of numerical evaluation discussed later in the paper. Furthermore, many problems can for a "first feel" 
be suitably approximated to be amenable to an analytic solution. Both of these features offer great help for intuition. 
At the same time, the Boltzmann equation can be derived from a full many-body quantum theory, for example with 
the Kadanoff-Baym equations as the starting point. This means that it in principle should be possible to keep track 
of all the approximation and their range of validity, and it also means that there is a way to work out more refined 



and extended equations, when needed. This stands in contrast to a derivation such as the one in section 111 However, 
whereas the Boltzmann approach is a very simple and instructive tool for simple problems such as deriving the Drude 
conductivity in an ordinary metal where spin is involved in a trivial way, it gets very complicated as soon as one 
considers more complicated systems such as graphene, where the pseudospin has a very non-trivial coupling with the 
momentum. Still, there is the possibility, as discussed in this paper, that a Boltzmann approach could make sense 
even for graphene-like systems and in some cases allow a fully analytical solution. This should be one clear motivation 
for investigating generalizations of Boltzmann equations in systems with non-trivial spin, in particular with some sort 
of spin-orbit coupling. 

In this paper it is shown that the simple but due to their chirality yet highly non-trivial model Hamiltonians (pi) 
with the generalization (pi) provide an ideal testing ground for scrutinizing and screening the kind of band-coherent 
generalizations of the Boltzmann equation that have been used in refs. IMTUl [2E1 [23 to discuss the conductivity in 
regime close to the Dirac point, that is, at the low charge carrier filling where the chemical potential is close to the 
zero of energy given by the touching of the two bands, the one with pseudo-spin up and the other with pseudo-spin 
down along the momentum. In this regime it might be relevant to not only consider the energy eigenstates of the free 
Hamiltonian, which describe the population of respective pseudo-spin band, but also consider electrons in quantum 
coherent superpositions of these eigenstates. (We will interchangeably use the terms spin-coherent and band-coherent 
to make it clear that we do not only consider populations of the eigenstates but also quantum superpositions of these.) 
The mentioned band-coherent Boltzmann equation approaches treat the spatial degrees of freedom semi-classically as 
in the standard Boltzmann treatment, but keep the treatment of the the band index — that is, of the (pseudo)spin — 
fully quantum mechanical. Because of the band-coherent terms, which are elusive to classical intuition, the collision 
integral of these equation cannot sensibly be given by an educated guess, but have to be derived from full quantum 
many-body theory. 

We are interested in to what extent these kinds of generalized Boltzmann approaches are in principal applicable, 
considering that the usual textbook criterion £kp 3> 1 for a Boltzmann treatment breaks down close to the Dirac 
point. We are also interested in if this approach is in practice possible, considering that for Dirac systems the collision 
integrals in these approaches usually become horrendously complicated, in particular if one takes the often neglected 
principal value terms seriously. The latter turn the static and uniform Boltzmann equation — usually a differential 
equation problem — into a integro-diffcrential problem. The mathematical challenges that come when solving these 
can be clearly seen in ref . |6]) . 

We show that the band-coherent Boltzmann equations considered in this paper are analytically solvable for the 
Hamiltonians (pi) and (p| with the impurity potential given by point-like impurities, and most remarkable the equations 
remain so with moderate additional complexity even when principal value terms are kept. The analytical form of 
these results makes the lack of equivalence of the different Boltzmann approaches conspicuous and facilitate a neat 
matching against other approaches that are considered more reliable close to the Dirac point, in our case a numerical 
evaluation of the finite-size Kubo formula. In the special case of bilayer graphene with point-like impurities, most of 
these approaches coincide. In this case the prediction of a band-coherent Boltzmann approach agrees well with the 
numerical evaluation, as shown in ref [26| . 

The source of this multiplicity among Boltzmann approaches comes from there being a wide variety of derivations 
of Boltzmann equations from quantum many-body theory, with both different starting points and different approxi- 
mations on the way. Even if the different approaches are sometimes motivated by different purposes, there are usually 
overlapping regimes where a comparison should be possible, like our case of weak and dilute impurities (implying that 
it should be enough to consider the lowest Born approximation and only non-crossed diagrams). In the context of 
graphene and spin physics the issue of their equivalence, or rather lack of such, was addressed by one of us [5]. It 
was shown that although there turns out to be agreement in some standard electronic systems, there is no general 
equivalence. In particular, for electronic systems where coupling between momentum and some spin type of variable 



is strong — as is the case for the here studies Hamiltonians with the pseudo-spin to momentum coupled term being 
the dominant term — at least six feasible and yet different band-coherent generalizations of the Boltzmann equations 
can be derived, and, as discussed in ref. [HI they all result in different spin-originated quantum corrections to the 
conductivity. On the other hand, in ordinary textbook systems with all degrees of freedom treated semi-classically, in 
particular for essentially spinlcss electrons with quadratic dispersion like in the ordinary metal, these derivations all 
agree, which is probably one of the reasons for the issue of equivalence of Boltzmann approaches not having received 
much attention. 

This abundance of candidates for the Boltzmann equation comes above all from the non-equilibrium Green's func- 
tions derivations. Here the Kadanoff-Baym (or Keldysh) equation is the starting point. In the derivation of a Boltz- 
mann equation one has to choose whether the real part of the self-energies should be related to relaxation (through 
the collision integral), as is the fate of the imaginary part, or if the real part instead goes into the renormalisation 
of the free drift of the semiclassical quasiparticle. The most important source of discrepancy, however, comes from 
the choice of Ansatz. The reason that one needs the discussed Ansatz is that one needs some approximation for 
expressing a two-time non-equilibrium Green's function G{t\ 1 t2) (other variables are here left implicit) in terms of 
of the one-time distribution function /(£) in order to get a closed equation for latter. The necessity of the Ansatz is 
illustrated in appendix [O in schematic explanation that has been abstracted and simplified from the kind of elements 
that occur in actual quantum many-body derivations of Boltzmann equations, see for example ref. [8] and references 
therein. 

The Green's function derivation can be contrasted with the more common and much more accessible quantum 
Liouville equation -based derivations of master equations (in our case the generalized Boltzmann equations) for the 
density matrix p(t) ~ /(£). Appendix |B] shows the origin of the Boltzmann collision integral in such a derivation. Here 
one is not confronted with the issues mentioned above. In the graphene related literature, almost all treatments have 
been of the Liouville type, see refs. IH1 13 HI 123 123 the to us known exceptions being refs. [HI EH In ref. [HI we showed 
that none of the existing Ansatzes in the Green's function literature gave a collision integral in the band-coherent 
Boltzmann equation that matched the collision integral we derived by a Liouville equation approach, and we invented 
yet another Ansatz to provide the missing link. With this we yet increased the number of Boltzmann equations that 
can be derived with a Green's function approach. 

One might think that a Liouville equation -based derivation is then the correct derivation as one does not seem to 
be confronted with these elements of choice. Our study does not support this view. Rather, we find that if any of the 
approaches can be ruled out, then it is the Liouville equation -based one, at least in its standard form. 

A central technical issue in the present paper is the role of principal value terms in the collision integral. They arise as 
byproduct in all derivations of Boltzmann equations mentioned above, also in the Liouville type of derivation. Principal 
value terms have been discussed in other contexts (under various names), for example the Boltzmann equations with 
spinless electrons subject to strong interactions, but have in the context of Dirac electrons and spin-orbit coupled 
electrons simply been left out, with the to us known exceptions being refs. £>] [HI The usual Boltzmann collision 
integrals have as typical ingredients a Dirac delta function 5{e s k — e|,) expressing the classically intuitive conservation 
of energy in a collision process. The quantum many-body derivations that lead to these expected collision terms will 

usually also result in similar terms, which instead of the delta functions involve principal values V f (e| — e|/) _1 ) . 

(Appendix [B| shows explicitly where they come about). The latter are rather elusive to a semiclassical interpretation. 
They might be connected with the Zitterbewegung property of Dirac electrons and spin-orbit coupled electrons. We 
believe that is motivated to address their presence and role seriously, in particular when the distribution describes 
quantum coherent superpositions of eigenstates. 

As mentioned above, the principal value terms make in most cases the Boltzmann equation horrendously complicated 
to solve. Maybe one reasons — on top of them being unintuitive — for them usually being neglected, most often without 
discussion. The special attractive feature of the models (2l) and ffih together with point-like impurities is that there are 
non- vanishing principal value terms due to the band-coherence but an analytical solution of the different Boltzmann 
equations is nonetheless still possible to infinite order in (Ikp) -1 in spin-originated quantum corrections and sums up 
to a neat result that remains finite even at the Dirac point ik-p = 0. Remarkable here is that the presence of principal 
value terms only moderately but not qualitatively increases the complexity of the solution. Furthermore, it turns out 
that the principal value terms modify the conductivity in a very non-trivial way. The studied setting is therefore a 
very interesting testing ground for screening of the different approaches and an addressing of the role and influence 
of principal value terms as well as the applicability of band-coherent generalizations of the Boltzmann equation. 



V. THE BAND-COHERENT BOLTZMANN TREATMENT 

In this section we are going to derive within the Boltzmann approach the pseudospin originated quantum corrections 
in the considered graphene type systems and in particular with principal value terms included. However, we have 
refrained from making the paper self-contained in this part, as this would make the paper far too ominous. For a 
technical understanding of the results familiarity is assumed with the derivation of Boltzmann equations in general 
as well as with ref. where technical details for the band-coherent case have been worked out in great detail and 
generality. This section functions as a guide for those readers who want to see how that work can be applied to the 
present context. Other reads might want to jump to the end of the section were the main results and conclusions 
have been gathered. 

In the Boltzmann approach the particles are described in terms of a distribution function f(x, k, t) living in classical 
phase space. In the presence of spin this function obtains a spin index, usually without any semiclassical approxima- 
tion, making / a density matrix in the spin indices. The Boltzmann equation then reads 

i[H, f] + d t f + i{v,, d x J} + eE t d k J = J[f] , (23) 

with the bracket { , } standing for anti-commutation. The presence of spin leads among other to the presence of 
the spin-precession term i[H,f]. However, spin-precession has also a classical analogue in the precession of angular 
momentum. The truly quantum mechanical ingredients are found in the interband components of {v,-,9 x< /} and in 
the collision integrals J\f] : for example in Pauli blocking factors and in the collision rates. 

The Boltzmann equation can be derived from quantum many-body theory, with the main issue being the derivation 
of the collision integral. One common derivation is to use the quantum Liouville equation as a starting point. A 
derivation of a collision integral can be found in appendix [B] Another starting point is the Kadanoff-Baym equations, 
from which one can derive 

[id t -H,G<} = E R G< - G<E A + E<G A - G R E< . (24) 

After Wigner transformation followed by gradient expansion and integration over the energy u> the left hand side 



turns into the left hand side of ( 23 ) , where by definition 

f(k,x,t) = f~G<(u,k,x,t), (25) 



An equivalent form of (251 in a non- Wigner transformed coordinates is 

f(xi,x 2l t) = G < (xi,ti,X2,t 2 )\t 1 =t 2 =t , (26) 

that is, the distribution function is related to the time-diagonal of the lesser Green's function. The right hand side of 



eq. (24) — involving self-energy terms S[G] which are functionals of the non-equilibrium Green's function G — turns into 
something that we can interpret as a collision integral. To derive a collision integral in the lowest Born approximation, 
the dressed Green's functions G R and G A are replaced by their non-interacting analogues G 0R and G 0A and in the 
self-energies the series in orders of the interaction is only kept to second order in the interaction and first order in the 
impurity density. The retarder non-interacting Green's function is in our case given by 



S=± K S—± 



with the projector 



S~ bs ■= ~(1 + »•*&*)• (28) 



The advanced Green's function is obtained by hermitian conjugation. 

As mentioned above, the collision integral is the ingredient of the Boltzmann equation that is the most demanding 
to derive. One reason is that it is non-linear in the distribution functions. It is not surprisingly in the collision 
integral that we will find differences between the different derivations. We will let the collision integral be divided 
up as J = J & + J v , with J 5 [/] containing the terms that contain delta functions and J7 P [/] containing the terms 
that contain principal value terms. In the case of monolayer graphene (Hamiltonian (pi) with N = 1) both parts are 
sensitive to the approach, even for point-like impurities. (See sec. 6 in ref. HI in particular eqs. (65,66).) For the 
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multilayer Hamiltonians with generalizations studied in the present paper, only the the principal value part will differ 
between the different approaches. 

Ref. [5] uses a nomenclature for the different approaches composed of a prefix and a suffix. The prefix Gl or G2 



indicates whether the collision integral derives from all self-energy terms as suggested in equation (24) or only from 



the imaginary part of the self-energies, namely from the right hand side of (the with (24) equivalent expression) 



[idt-H -Re £ R G<] - [£<, Re G R ] = i{ImE R ,G < } - i{Y,<,ImG R } . (29) 

In the latter case the real part of the self-energies to the left are instead assumed to renormalize the free drift. (As 
ref. [5] was not able to carry out this renormalization in the spinful case, the statements and results for the approach 
G2 have to be taken with more caution than those of Gl.) The suffix GKBA (Generalized Kadanoff-Baym Ansatz 
[28]). A A ("Alternative Ansatz" or " Anti-ordered Ansatz" |8J and SKBA (Symmetrized Kadanoff-Baym Ansatz [8]) 
refers to which Ansatz has been used. In the Green's function derivations considered in this paper the issue of Ansatz 
is the central one. 

The Ansatz gives an explicit relation between the non-equilibrium Green's function G and the distribution function 



/. By definition one has f(t) = G < (t, t) , or in the Wigner-transformed language equation (25). This implicit relation 



is enough for the left hand side of (24) to turn all Green's functions into distribution functions and this thanks to 
the left hand side being linear in Green's functions. However, the right hand side is non-linear in Green's functions, 
wherefore the implicit relation such as ( 25 ) is not enough, as explained in appendix O Here one needs an explicit 



relation between / and G < (and one needs explicit expressions for G 0R and G 0A , here given by (27)). Such a relation 
between a semiclassical object and a fully quantum mechanical object has of course to be an approximation. The 
Ansatz is the approximation that provides the explicit relation between / and G < , allowing one to translate also in 
the other direction, that is, finding G < if one knows /. This is necessary for distilling a closed equation for / — as is 
the Boltzmann equation — out of equation (24). (See again appendix |C| for a more detailed explanation. 



There are several possible candidates for the required approximation involved in the Ansatz. The considered 
Ansatzes can for G < (xi,ti,X2,t2) in the lowest Born approximation be written as 

GKBA G<(h,t 2 ) = J G 0R (t 1 ,t 2 )/(t 2 )- l /(i 1 )G 0A (t 1 ,i 2 ) 

AA G<(h,t 2 ) = l /(ii)G 0R (t 1 ,t 2 )-*G 0A (i 1 ,i 2 )/(i 2 ) (30) 

SKBA G<{hM) = |A°(ti,t2)/(t2) + |/(ti)A (ti,t 2 ) 

where A = i(G R — G A ) is the spectral function. All Ansatzes are consistent with /(£) = G < (t,i), but go beyond 
the definition by providing information also about the timc-off-diagonal part of G < . Note that the Ansatz SKBA is 



the symmetrization of GKBA and A A. In equation (30) the spacial two-point indices (a;i,x 2 ) have been suppressed 
for convenience, and so has on the right hand side also the convolution product in these. (In the time indices there 
is no convolution product.) For the background to these Ansatzes, as well as for the approximations that follow 
upon the Wigner transformation and with the Born approximation, see ref. HI Here we just remention that the 
combination GlAA is what provides the missing link between the Green's function derived collision integral and the 
Liouville equation derived one (as derived in appendix [B]), at least in the lowest Born approximation, which is what 
we consider. Let us also stress that the list of approaches need not be exhaustive, but provides only a set of natural 
candidates. 

The principal value part J v [/] of the collision integral can be written as J v — J FX + J VY where J px [/] is the 
part that does not change sign when one compares the approach G1GKBA to the approaches GlAA, whereas J VY \f\ 
is the part the changes sign. The principal value terms J v for the other approaches can be composed from these two 
parts and is given by the table 

X Y 



Gl GKBA + + 

Gl AA H 

Gl SKBA + (31) 

G2 GKBA + 
G2 AA - 
G2 SKBA 0. 

Using the matrix decomposition in eq. QSp the Boltzmann equation including the collision integral will in a linear 
approximation in the electric field decompose in the Fourier modes of the angular variables 9 of the momentum k. 
For the terms J p [f^ E '], with f( E > being searched for linear-in-electric-field correction to the known equilibrium part 
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/ cq one has according to ref. [5] (for notational convenience the superscript (E) is suppressed; the prime superscript 
indicates that the component of / is a function of k' and not of fc)|33j 



with 



qPX 
rPY 



bl 
PX 



•-'ex 

qPY 

rrPX 
<-*zl 

q-PY 

<->zl 



= 



q-PY _ 



+i k> ^if v - sin N c A # sin A ®fzi , 

-Ik' "%T (-V^cosN c A0f zl +V + cosAdf zl ) 

-Ik' zf (-■P+f z i + V-cosN c AecosAef zl ) 

-Ik' %^ {V-CosNcAefa-V+iismNzAdsinAdf^ + cosNcAdcosAefa)} 



(32) 



V± 



± 



6 + 6' 6-6' 



(33) 



In the case of monolayer graphene, both J px and J VY are typically non-zero, and the six candidates in (311 lead 
to six results — all different — for the leading quantum correction. In contrast, for the case of N c > 2 (including the 
multilayer Hamiltonians Q with N > 2) combined with point-like impurities iW kk ' — 2ttu^ =: Wq) most of the terms 



in (32) vanish trivially, namely all the terms that contain trigonometric functions. This leaves as only nonzero terms 



w t 



qPY „ ,jr f vv kh> r> 

J ex - +JziJk' ~^T' + -■ 

n-PY ._ r f w kk' -n 



2tt 



+/ 2 iA fc 
— /ciAfc . 



The vanishing of all terms J implies that we cannot distinguish between the approaches G\ and G2. We henceforth 
assume that we deal with Gl and suppress this prefix. More interestingly, the fact that some terms J VY survive 
means that we can distinguish between the different Ansatzes GKBA, AA and SKBA. Furthermore, the integral in 
Afc can be worked out analytically as (see appendix |A| 



Ak = 



2W 

(2tt) 2 6 



k'dk' 



V 



cot 



N d 



Ttr 



2-N d 



(34) 



Note that the derivation requires Ad > 1. 

Since the surviving principal value terms only involve fk but not fk> we can plug out the distribution function out 
of the integrals, which allows us to turn the integro-diffcrential equations into differential equations, as is typically 
the case for static and uniform Boltzmann equations. Thanks to this we essentially have the same type of Boltzmann 
equation as when the principal value terms were neglected. Section 6 in ref. |8] treated the latter problem in a general 
setting and the results can be easily translated to the present case. Equation (69) in ref. |8] is in the present case 
replaced by the following equation 



E„ — iE„ 



( 4/o q \ 

1 k Jf, 

V o / 



( 1 
1 







1 







26-j^A 



AE) 

J bX 



\ i/A - 26 1 J 



AE) 

Jcl 

V /if } / 



(35) 



Most importantly one has the simple shift 26 — > 26 — vh in the strength of the precession term due to the principal 
value terms, with the sign depending on Ansatz according to 

GKBA v = +1 
AA v = -1 
SKBA v = 0. 

Furthermore, N — > N c in the right hand side. The chiral origin of that N can be found above equation (65) in ref. |8] 
Other changes are that I + , I x and 1 K coalesce to X + , in the present paper denoted I, and that X s = in the present 
setup. The result in equation (75) in ref . [5] becomes in the present case 



a 



CD 

j-N c 



J£d* 



b 

k (2b- 



(36) 



i/A) 2 +I 2 
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The result that there are no quantum corrections to a 1 gives an explanation from the perspective of the Boltzmann 

made sense. There it was assumed that <7 lnter 



III 



approach to why one of the central assumptions in section 

all the quantum corrections due to interband effects, whereas cr lntra is identical 



to 



contains 
the Drude conductivity ud, so 



contains no corrections. The contribution cr lntor relies on the term v lnter . However, according to definition of current 
in terms of the Boltzmann distribution fucntion in equation (M) , the velocity component v mtcr is what the distribution 
function component /g combines with. In the present case this component authors singel-handedly the correction 
a 11 . To see this, look at eq. (72) and eq. (75) in ref. [5] and remember that I s = in our case. Note that for the 
monolayer problem there is no such clean separation. There also /g, which contributes to the current through v lntra , 
contains quantum corrections, wherefore we expect that for the monolayer case the crucial assumption in section III 
not to be valid. 

Even in the presence of principal value terms the correction term a 11 can be rendered on a neat analytical form 
With e and h reintroduced, we arrive at the main result (see appendix [A]) 



II 



e 2 £k F 


e 2 N d b F T tr 




ffD - h 2 " 


h 2h 




e 2 N c 2 


arctan 

\2 \ 


'2lk F 


h 8(7V d - 1) 


v N d 



J/ cot ■ 



N r] 



(37) 



For v = and N c — N d this is the result in equation (78) in ref. [51 
The conductivity at the Dirac point is found by setting £k F = 0. 
since a 1 = 0. Using that 



The conductivity is then only composed of a 



ii 



arctan cot 



( cot w) 



N-2^ 

2N '' 



(38) 



one has 



a(£k F = 0) = 



h 2 



\N d + v{N d - 2)) 
8(7V d _ i)jV d 



h 8 



Nd 2 

2(N d -l) 

N c 2 
N d (N d -l 



GKBA 
SKBA 
AA. 



(39) 



Here we gather the main points of this section 



That all terms J vx in equation (32 1 are zero means according to the scheme (31) that the multilayer problem 



with point-like impurities docs not allow to distinguish between all the approaches referred to in table (31) . 
However, interestingly it is J VY that survives, which according to (31) allows us to distinguish between the 
three Ansatzes GKBA, AA and SKBA. 

• In contrast to the monolayer problem, the surviving principal terms depend only on fy, but not on /&/. One can 
therefore plug out the distribution function out of the principal value integral, leaving a momentum dependent 
integral A& that is independent of the distribution function. We are therefore luckily not forced to solve the kind 
of integro-diffcrcntial equation considered by Auslender et al. [B] , but only the differential equations of the level 
of complexity that we have without principal terms. This is an appealing feature of the multilayer problem (and 
the general Hamiltonian (p| with N c , N d > 2) combined with point-like impurities. (Note that if the impurities 
are not point-like, then the terms with trigonometric functions in (32) normally do not vanish and we are stuck 
with the full integro-differential equations.) 

• An interesting special case is the one with quadratic dispersion N d — 2 but general chirality N c > 2. Then A = 
and the principal value terms vanish all together, with all approaches agreeing, see figure [TJ This is the fortunate 
circumstance that allowed us to neglect principal value terms and the question of different approaches in the case 
of bilayer graphene with point-like impurities studied in ref. 1261 It also results in density of states independent 
of filling, which we find is a big advantage for the numerical analysis of the Kubo formula as performed in ref. 1261 
and in the next section of the present paper. One important reason for introducing the artificial Hamiltonian 
^ is to be able to study the role of increasing chiral winding numbers in such a special setting. 

• The Boltzmann derivation supports the crucial but loose assumption made in section [TTT] that all the quantum 
corrections are related to the interband velocity component v mtcr , contributing to er mtGr , whereas cr lntra , paired 
with the velocity component v lntra , comes with no corrections. See the discussion under equation (36). 
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It is interesting to note that at £kp = the conductivity derived with the Ansatz GKBA in the Boltzmann 



approach matches the conductivity ( 20 1 derived in section. Ill with the heuristic analytic treatment of the Kubo 
formula. However, for £kp > the two results no longer agree except in the special case of quadratic dispersion 
Ad = 2, see figures [T] and [2j Thus, in general the heuristic derivation finds no equivalent among the discussed 
Boltzmann derivations. 



• As mentioned already in section. Ill the conductivity a(£kp = 0) is also the conductivity minimum when 
N c < Ad. However, it is not necessarily so when A c > Ad, in which case the conductivity minimum can occur 
at some strictly positive £kp whereas £kp corresponds to a local maximum in the conductivity as a function of 
£kp, see figure [TJ 

• It is interesting that the result derived with the Ansatz AA, which would also be the result derived with a Liouville 
equation derivation, predicts a conductivity minimum that decreases with the number of layers A(= A c = ./Yd) 

2 

of rhombohedrally stacked graphcnc, reaching cr m i n = xf (half of the bilayer result) in the limit A — > oo. The 
other Ansatzes and the heuristic Kubo derivation predicts that the conductivity minimum increases with the 
number of layers, see figure [2j 

As one further point we mention that based on particle-hole symmetry we would expect the curve to be mirror 
symmetric if we were to plot a as a function of charge carrier density n (which is negative for holes). However, since 
for positive n we have n ex £kp in the studied setups, this would mean that if the conductivity at a(£kp) is not 
horizontal at £kp = 0, then a{n) has a kink there where it passes n = 0. However, such a kink does not seem very 
feasible and the numerical analysis does not suggest any kinks. As for this issue the analytical Boltzmann treatments 
discussed by us should not be expected to give a valid description at £kp = 0. It is interesting to note that the 
heuristically derived result based on the Kubo formula produces no such kink for N c = Ad. 

VI. NUMERICAL STUDY OF THE KUBO FORMULA 
A. General issues 



In this section we present the results obtained by a numerical evaluation of the Kubo formula ( 12 1. The Kubo formula 
can be applied to a finite system connected to leads by choosing r\ — g^Se. The finite r\ expresses the broadening of 
levels in a finite system given by a possibility of the particle to escape from the system and is related to the average 
energy spacing 5e = (DL 2 )^ 1 at the Fermi energy. Here, D is the density of states at the fermi surface. In the 
non-interacting case it is given by 

D = W/ 1 ""' (40) 

i?t is the dimensionless Thouless conductivity. In principle (?t(£f) can be estimated. However, in the case that system 
is dirty enough so that the collision time To, which in our case equals r tr , is much smaller than the Thouless escape 
time 7"Th = ^/(9T^ e )i it turns out that, whereas the conductivity goes to zero for g^ too small or too large, it is 
rather insensitive to the precise value of <?t in an broad intermediate range of the latter. As we are interested in 
bulk conductivity and do not want to probe the finiteness of the system and the presence of leads, it is exactly the 
conductivity in this t/x-insensitive plateau that we are after. As in ref. |26l I29| we typically set <?t to be about 10. 
The condition for a dirty system can also be written as 

T Th fiDL 2 2tt Dn % u 2 ■nD 2 N i u 2 , , 

i« — = ir^— = — ( 41 ) 

Ttr 5t n 2 g T 

i.e. 



We see that to assure ( 42 1 , the impurities should be strong enough, or the number of them has to be big enough. 
On the other hand, since we are focused on comparing with Boltzmann type equations that were derived in the 
lowest Born approximation, it is important that the interaction strength is chosen small enough so that higher Born 
corrections can be expected to be negligible. Since higher Born corrections typically come with factors of 
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the leading contributions will for point-like impurities U kk > = uq/L 2 come from the smallest denominators e k — e k > . 
An estimate of the average minimal size of this difference is given by average energy spacing (DL 2 )^ 1 . We have 
therefore the second condition 

u \ 1 

U kk > = J2 « Se = ^2 =*• U « J) ( 44 ) 



Putting ( 42 ) and ( 44 1 together we obtain the following double condition on the parameters 

5t 



, r <uL>«l. (45) 



We see that we need to have y/gT/irNi -C 1, which is rather easy to satisfy. The challenge is to keep uD within the 



two inequalities in the condition ( 45 ) . For ordinary point-like impurities the meeting of this condition will turn out 
to require an individual adjustment for each different N&. For the sake of comparability of different N we have, for 
reason soon to be explained, ad hoc introduced an artificial modification of the point-like potential, with the strength 
being dependent of kp. The two point-like potentials u which we will study numerically are 

u _ f JVd7 2 7r dL c t ordinary potential ,^, 

1 H artificial potential . 

Here Ci is a dimensionless constant and k — 2nh/L. For N^ = 2 the two potentials coincide. 

For a general dispersion the density of states D depends on kp. In the case of bilayers (N = 2) studied in ref. [25] 
the situation was particularly simple, because for a quadratic dispersion e ex k 2 the density of states D is a constant. 



A set of parameters satisfying the condition (45 1 for some kp satisfied the condition then for all kp. In contrast, in the 
multilayer case with the dispersion |e| oc fc*^, the density of states D depends on kp. For the ordinary potential with 
a constant strength this implies, that for low or high enough kp left or right part of the condition will be increasingly 
violated. At the low enough kp the assumption of no higher Born corrections is violated, at the high enough kp the 
assumption r tr <C rxh of the system being dirty is violated. The strength of the ordinary potential has therefore to 
be adjusted with a great care with respect to the power of dispersion as well to interval of kp that we wish to study. 
We have introduced the artificial potential to carry over the fortunate circumstances for bilayers to all N: with the 
artificial potential the condition ( [45] ) is satisfied for all kp equally and this for all N such that the other parameters, 
for example Cj, do not need to be adjusted as we increase N. 

We define the dimensionless momentum q — (hj n)k with its special value N e = irqp. The impurity potential then 
reads 



u _ Nd'JK d J a ordinary potential 

L 2 2tt 1 aq F d ~ artificial potential 



U k = T2= ^ <i . JV d -2 _„^„.. 1 _. i .^_ 1 ' ( 47 ) 



The functional dependence of £kp on the parameters for the two types of potential is 

_ / c?m q F Nd ~ 2 ordinary potential 

I -r^rQv = -t-to* = -^i— artificial potential 

Note that with the artificial potential we recover the bilayer situation that £kp ex n/rii, which is not the case for 
ordinary point-like impurities in the multilayer case. 

Independent on the setup of the potential strength there is an additional upper constraint on (.kp coming with 
the momentum cutoff introduced by the necessarily finite matrices in the numerical treatment. The above defined 
dimensionless momentum q is for a finite system integer valued and label the momentum states pi = 2j ^ k qi- For 
matrices of rank 1q\ (2 coming from the spin) we have to keep ikp < l\k\, where the cutoff is given by 

{c-f Q AT n 

-^f^x d ordinary potential 
-^jrq\ artificial potential 

B. Numerical results 

We will now discuss our numerical results presented in Figures [3] and HI The plots show the conductivity in units of 
e 2 /h as a function of the dimensionless variable (.kp. We consider two cases of disorder: ordinary point-like impurities 
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N = 3 



N =4 



N = 6 




FIG. 3: The numerical finite-size Kubo conductivity computed for the chiral particles for a given N(— Nd — N c ) in the 
presence of ordinary point-like impurities (red squares) and an artificial scattering potential (green triangles), respectively. 
The solid curves show the results for the different analytical Boltzmann approaches, distinguished by the Ansatz used for the 
derivation of the collision integral. The analytical curve in the middle corresponds to the Ansatz SKBA in which principle 
value terms are absent in the collision integral. Of the compared Boltzmann approaches this one appears come closest to the 
numerical results in the range of low momenta. The parameters used for the conductivity computations in the presence of 
the artificial potential are q^ = 14, a — 0.4, #t = 12 at all N^. With the artificial potential the obeying of condition (451 is 
independent of the value of fci? as well as of N. In the case of the ordinary point-like impurities, the parameters must be chosen 
different for each N^ to satisfy condition (J45|. We have used JVd = 2: q\ — 14, a = 0.4, g^ — 12; iVd = 3: q\ — 12, Cj = 2.4, 
ff T = 12; iV d = 4: g A = 12, c, = 10, g T = 12; iV d = 6: qx = 10, c, = 200, g T = 12. 



N = 2 



N =4 




FIG. 4: The numerical Kubo conductivity (red squares) computed for the chiral particles with N — 2 and N = 4, respectively, 
in the presence of ordinary point-like impurities. Also shown is the intraband component a 1 a (green circles) and the interband 
component u mtor (blue triangles) The qualitative behavior is matched by the analytical formulae derived within the Boltzmann 
approach with the Ansatz SKBA. The parameters used for the conductivity computations are the same as in the previous 
figure. 
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as well as the artificial point-like impurities, as introduced in Sec. |VI A| Figure [3] shows the total conductivity for 
N = AT(j = N c equal to 2, 3, 4 and 6, with the corresponding analytical dependence given by the solid curves. The 
numerical results are consistent with the analytical prediction that the conductivity does not fall down to zero at 
£kp = 0, i. e. the main qualitative feature - the finite conductivity minimum at the neutrality point - is captured. 
Figure H shows the intra- and interband contributions cr lntra and <7 mtcr . The intraband contribution tr lntra (£fcF) is 
monotonically increasing function whereas the interband contribution <r mter (^fcF) is a monotonically decaying function. 
In our study we checked that the values a mtr& + a mtcr are in most cases indistinguishable from the total conductivity. 
(For N = 2 we actually see a visible but still small discrepancy.) This indicates that there are no other relevant 
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contributions to a than cr ,ntra and o- mtcr as defined in Eq. 

From the theoretical treatments of the present problem we expect er lntra (^£;F) to contain no interband contributions. 
Therefore its value should be close to the Drude conductivity od = (e 2 /h)£k F /2, at least for large enough £k F , unless 
there are deviations due to for example weak (anti-)localization. For N = 2 the interband contribution clearly follows 
the Drude prediction, as seen in the left panel of Fig |4] Initially it surprised us that for higher N there is a clear 
discrepancy between the numerical results and the Drude prediction, as can for example be seen in the right panel of 
Fig [4] Subsequently we have understood the reason for this and hope to discuss it in detail in upcoming work. Like 
any semiclassical approach the Boltzmann approach, and in particular the derivation of the Drude conductivity, has 
as a prerequisite a unique, well defined quasiparticle peak. For higher TV this condition is approximately met at low 
£k F , but no longer at higher £k F . Interestingly, this means that if the Boltzmann approach can at all be applied to 
multilayer hamiltonians, then only close to the Dirac point. 

The results in Fig. [3] show that although the artificial potential is rather different from the ordinary potential of 
point-like impurities, the two potentials give very similar results. The discrepancy between them for higher N at large 
enough £k F should have the same reasons as discussed in the previous paragraph. When the Drude picture fails, we 
no longer can expect our translation between fep and £k F to be correct. 

We now come to what the numerical results imply for the choice of Ansatz. In the case of N larger than 2, the 
analytical approaches disagree and show even qualitative differences. For example, the Boltzmann approach with the 
Ansatz AA (equivalent to the Boltzmann equation derived with a Liouville equation approach) predicts a conductivity 
minimum that decreases with N whereas an increase with N is observed in the other considered approaches. When 
comparing the numerical results obtained from the Kubo formula with the analytical approaches, we realize that the 
agreement is mainly on the qualitative level. However, the qualitative features appear to be distinct enough for us to 
be able to rule out approaches. Fig. [3] reveals the numerical and analytical conductivity curves for different values of 
N. As the numerical results clearly point at the conductivity minimum increasing with increasing N we would rule 
out the Boltzmann approach with the Ansatz AA and have therefore left it out in Fig. [4] not to overload the figure. 
As the numerical results lie closer to the prediction provided by the Boltzmann conductivity with the Ansatz SKBA, 
we only include these and also leave out the results from the Boltzmann conductivity with the Ansatz GKBA. 

For our methodological studies, we have also considered the problem of fixing the power of the dispersion in 
Hamiltonian p| to be quadratic (JV<j = 2) and only varying the chiral winding number N c . This artificial model 
helps us to distinguish between how the chiral winding number influences the conductivity from how the power of 
dispersion does so. The latter factor can also be thought of as the influence of the density states, as increasing 
N^ leads to an increased density of states at the Dirac point (whereas far away from the Dirac point there is the 
opposite effect). Furthermore, since all the by us considered analytical approaches agree for a quadratic dispersion 
this setting is particularly good for obtaining an impression of what kind of qualitative and quantitative agreement we 
in general can expect between the analytical and numerical approaches. As seen in Fig[5j the qualitative agreement 
is good. In particular, the local maximum of conductivity at £k F = increases with N c in agreement with the 
theoretical prediction in Fig [l] However, the quantitative agreement close to the Dirac point becomes increasingly 
unsatisfactory with increasing N c . The rightmost panel in Fig [5] shows that the quantitative discrepancy between 
the numerical and analytical results observed in the total conductivity is found in the interband component. For the 
intraband component, on the other hand, the quantitative agreement is very convincing. In particular, we do not 
find any dependence on N c in the intraband conductivity. We mention that we observe the interband component to 
a convincing degree to obey a N c 2 scaling. 

For both considered problems we can conclude that the analytical approaches agree qualitatively with the behavior 
in the numerical results, but quantitatively the agreement is still unsatisfactory close to the Dirac point. 

VII. CONCLUSION 

In this paper we have studied the conductivity of chiral particles within the two-band model suitable for the description 
of low-energy electron excitations in idealized rhombohcdrally stacked multilayer graphene. Our aim has not been to 
capture properties of a realistic graphene sample, but to use a clean, idealized setting for calculating the conductivity 
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FIG. 5: The total Kubo conductivity (left panel) and its components cr mtra (middle panel) and <j lntor (right panel) computed 
for the chiral hamiltonian with JVd = 2 but with chiral winding number N c varied. The solid curves represent the corresponding 
analytical results. Note that when JVd = 2, the principle value terms do not play a role and all approaches agree, as shown in 
Figfll Furthermore, for N^ = 2 the artificial potential and the ordinary potential are identical. The results for different iV c are 
easy to compare on the same footing since the validity of condition (45 I is independent of /cf and N c . The parameters used are 
Cl = 0.4, 3T = 12, q x = 12. 



in the presence of chiral interband effects, and use this idealized setting for comparing analytical approaches against 
each other. To compute the conductivity we have employed the finite-size Kubo formula and performed an exact 
diagonalization of the total Hamiltonian consisting of the mentioned chiral kinetic Hamiltonian and a point-like 
impurity potential. Among the analytical approaches we both presented a heuristic derivation based on the Kubo 
formula as well as studied a pseudospin-coherent generalized Boltzmann equation. In the latter we stressed the leading 
role of the chiral winding number N c in the band-coherent conductivity and showed that the choice of Ansatz has a 
huge impact on the results close to the Dirac point, see Figures [T] [2j 

In the numerical analysis we saw that the conductivity minimum increases with the layer number N, see Figure [3] 
Comparison with the analytical results rules out the band-coherent Boltzmann equation derived with a Green's 
function approach using the Ansatz AA (the one with v = — 1 in the principal value collision integral terms, see 
Figure^]) which predicts the opposite behavior. The band-coherent Boltzmann equation is equivalent to the one derived 
within the Liouville equation approach, which is the simplest and most popular approach of deriving generalized 
Boltzmann equations in the spin-orbit interaction context. The conclusion made above therefore rules out the Liouville 
equation approach as well. This approach seems to be unsuitable in general when band-coherent effects become 
important, although it might be acceptable as long as the principle value collision integral terms do not play a role. 
The evidence against the Liouville equation approach for the derivation of the principal value terms in the collision 
integral is one of the main and most surpising results of the present paper, especially if one considers the simplicity 
and apparent self-sufficiency of the Liouville-like derivation. 

On the constructive side our hope was to single out one approach for deriving band-coherent Boltzmann equations. 
The Green's function derivation with the Ansatzes GKBA [y = 1) and SKBA (y = 0) result in analytical predictions 
closer to the numerical results than the mentioned derivation using the Liouville equation. They predict correctly 
a conductivity minimum which increases with the number of layers, see Fig [3] The best results have been obtained 
within the approach using the Ansatz SKBA. In the corresponding collision integral the principal value terms are 
absent by exact cancellations. The discarding of the competitor GKBA would also entail the discarding of the 
Liouville-type derivation with the Markov approximation done not in the interaction picture but in the Schrodinger 
picture, see Appendix [B] Support for the Ansatz SKBA suggests that principal value terms are not a part of the 
correct band-coherent generalization of the Boltzmann collision integral. In an approach were they are present one 
does best in neglecting them. One example of this is the appendix in ref. 1301 were a Green's function approach with 
the GKBA is used. The principal value terms found in the derivation are neglected, without detriment to the end 
result. 

One interesting observation was made when we decoupled the chiral winding number N c from the power N^ of the 
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dispersion, and kept the latter fixed at N^ = 2, in which case all analytical approaches agree. As seen in the leftmost 
panel of Figure [5] the total conductivity shoots up when £k F approaches the neutrality point. Moreover, the value at 
£kp increase with the winding number N c , as predicted by our analytical expression. For a given N c the value at the 
Dirac point is larger than in the case when TVa = iV c , shown in Figure [3j This supports the analytical prediction that 
an increase of the winding number N c increases the conductivity close to the Dirac point, whereas the increase of the 
power of dispersion N<± has the opposite effect. The latter effect is somewhat surprising considering that increasing 
iVd leads to an increased density of states at the Dirac point. In the multilayer case N c = N^ = N these two effects 
partially cancel each other, with the net effect still being an increase of the conductivity at the Dirac point with 
increasing N. 

At higher values of £fcp, the numerically calculated conductivity for higher layer number N showed a feature that 
we did not expect from the analytical treatment. The total conductivity is found to be significantly below the Drudc 
conductivity at higher electron densities, although one would imagine that one is safely in the Boltzmann regime 
Ik-p 3> 1 far away from the Dirac point. We have found that this discrepancy between the numerical Kubo and 
analytical Drude predictions is due to the lack of a unique quasiparticle peak and can be eliminated if the scattering 
potential has not a point-like structure. The problem is well known for the delta-functional scattering potential in the 
two-dimensional case. (See e.g. the review [31j.) The Born approximation tends to underestimate the scattering cross 
section in this case. Since this issue does not seem to be related to interband coherence and to chirality but rather to 
the properties of spectral functions of excitations in a given system, we postpone the discussion to subsequent work 
dealing with the electron transport of chiral particles far from the neutrality point. 

To conclude, we have shown that the band-coherent Boltzmann approach has a potential of being a theory of 
conductivity of chiral particles close to the neutrality point, even though this approach seems to be formally invalid 
in this region. For the derivation of Boltzmann-like band-coherent kinetic equations the Liouville equation approach 
should be used with great precaution. In particular, it turns out to be unsuitable for the derivation of the principal 
value terms in the collision integral. However, it can probably be used for the qualitative analysis of the band-coherent 
conductivity as long as these terms play no role. 
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Appendix A: Calculations used for the band-coherent Boltzmann treatment 
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Appendix B: Derivation of the collision integral from the quantum Liouville equation 

The simplest derivation of a collision integral is probably the one in which the quantum Liouville equation is iterated 
to second order in the interaction. Let H to t = Hq + V where V is an interaction switched on at a time to in the 
remote past. The quantum Liouville equation for the density matrix p SlS2 (t 1 x\,X2) in the interaction picture is 

id t p l = [V\t) lP l ] (Bl) 

with A 1 {t) = e lHot A(0)e- iH(,t =Ul(t)A(0)U (t). This is easily integrated to give 

p\t) = p\t,)-ifAt'[V\t') lP \t')] (B2) 



which, when inserted back into (Bl), yields 



d t p l = -i[V\t), P \t )}- / dt'[V\t),[V l (t'),p\t 



In 



i^y^]- f &t>[V\t),[V l {t>),p\t)]]~ l f f &t>[V\t)AV l {H)AV\t"),p l {t")]] 



(B3) 



after a first and second iteration, respectively (the second time with the interval of integration (to,t) replaced by 
(£',£) ). So far the equations are exact. The Born approximation allows us to get a closed equation for p at time t to 
second order in the interaction V . To this end we remove the last term in the second row. Equivalent to this is to 
replace in the last term of the first row, the full evolution with the free evolution, i.e. let p l (t') = p l (t), which has the 
appearance of a Markov approximation. Back in the Schrodinger picture the assumption of free evolution reads 

p(t') = U (t' \t) p{t)Ul{i! \t) = e - lHa{t ' V p{t) e lHo{t ' ' t] (B4) 

and after reorganizing evolution operators one obtains the kinetic equation in the Schrodinger picture, 

d t p{t) + i[H 0l p{t)] = ~i[V,e- lH ^~ t ^p{t )e lH ^ t - t ^]- [ °dT[V,le- tH ° T Ve lH ° T ,p(t)]]. (B5) 

Jo 

So far, this kinetic equation is locally time reversible (globally not, since we switched on the interaction). To capture 
the decoherence due to other processes (phonons etc.) we do not want the state to depend on correlations in the 
remote past. Therefore we include a factor e~ ? ' T in the integral to impose this loss of memory. This introduces time 
irreversibility. This factor also regularizes the integral and allows us to send t — > — oo. In translating the evolution 
operators into Green's functions, the last term can be written as 

J[p] = -J^[V,[G™VG° A ,p}}, (B6) 

where we anticipated that the last term will become the collision integral J . 

Another source of irreversibility comes with the impurity averaging procedure, a coarse graining that also captures 
the decoherence due to phonons, for example. Then V in the first term becomes just a number ~ V(r = 0) and 
the commutator vanishes. |34j For the terms linear in 7ii mp one finds for example (summation over repeated indices 
implicit) 

(VG 0R VG 0A p) kk , — > n imp S(k -ki + hi - k 2 )u kkl Giy klk2 G° k A Pk2k , = ^G° k A Pkk ,, (B7) 

where we introduced the retarded self-energy E^ = ni mp (uG 0R u) k . However, in a term like 

(V P G m VG 0A ) kk , — > n imp S(k -k 1 + k 2 - k')u kklPklk2 G™ 2 u k2k ,Gl A (B8) 

the delta function seems to offer no simplification at all. 

At this point we can attain further simplification if we say that in the collision integral we are not interested in any 
contributions which have to do with non-diagonality in momentum. This actually amounts to saying that we are not 
interested in any gradient expansion corrections to the collision integral: 

(VpG m VG 0A ) kk — > n- mp u kk , Pk ,G^u k , k Gl A . (B9) 
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The collision integral can then be written as 

J[p(k, x, t)} = - J W kk> J J^( Pk G k *G° k $ + G° k ?G k A Pk Pk ,GfGi A G° R G#V) (BIO) 



with the transition matrix W kk > = 27mi mp |u fcfc '| 2 for spinless impurities. The collision integral (BIO) is identical to 
the one derived with the Green's function approach GlAA. The collision integral will in general contain both delta 
function terms and principal value terms. 

Note that if one were to adopt a Markov approximation instead the Schrodinger pciture as p(t') — > p(t) (as done for 
example in ref. [SJ rather than p l (t') — > p l (t) as in the Dirac picture, then the evolution operators cancel each other 
out in a different way so that at the end they sit around the entire inner commutator rather than only around the 
inner V, 

J[p] = - [^[V,G™[V,p}G° A }. (Bll) 



This last Markov approximation is not equivalent to neglecting the 0(V 3 ) term in the last line of (B3), that is, it is 
not equivalent to what we understand as the Born approximation. What is a correct Markov approximation is not 
important to us. The base line is that we want to derive the collision integral in the lowest Born approximation, and 



for this the steps leading to equation (B6) seem to us to be the right ones. 

The different structure of of the double commutators in (B6) and (|B11) are crucial. The collision integral (Bll) 



coincides now with the one derived with the Green's function approach G1GKBA. 

Appendix C: Simple example of why one needs an Ansatz 

Suppose that by definition 

duj A(uj,t,x,p) = a(t,x,p) 

dujB(uj,t,x,p) = 1 (CI) 

for some quantum correlation functions A and B. Although there is a simple relation between A(uj, t, x,p) and 
a(£,x,p), the former will typically contain much more information than the latter, and some of the information goes 
lost when we integrate over the variable uj. 

Suppose there is an equation of motion such as 

d t A(w, t, x, p) + v p ■ d x A(oj, t, x, p) = . (C2) 

We can integrate the whole equation over uj to get rid of this variable and used the definitions to obtain 

d t a(t, x, p) + v p -d x a(t,x, p) = 0, (C3) 

that is, a closed equation for a(i,a;,p), meaning that the original unknown quantity A only enters in the form of a. 
The latter equation of motion contains of course less information than the former (some of the information was lost in 
the to- integration) , but might be that the information in the latter equation is just about the amount of information 
we are interested in and that the equation is now analytically more tractable. 
However, for an equation like 

d t A(u;,t,x,p) + v p -d a .A(w,t,x,p) = / Y pp -B(ui,t. x,p)A(uj,t,x,p') (C4) 

J P ' 

we have for the right hand side no use of the definitions when we want to integrate the whole equation over uj to get 
rid of this variable. The definitions tell us only what the outcome in terms of w independent quantities are when 
A(to,t, x,p) and B(ui,t,x,p) are standing alone or only together with quantities that do not depend on u. There is 
no general identity that would give us a result of the integral 

dui B(uj,t,x,p)A(uj,t 7 x 7 p). (C5) 
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To make progress we need more explicit information about A(tu,t,x,p) and B(ui,t,x,p) than what is contained in 
the definitions. For example, we might know that we can work with 

A(uj,t,x 7 p) = a(t,x,p)5(uj — e p ) 

B(uj,t,x,p) = 6(u-€p). (C6) 

These relations are consistent with the definitions but provide us with more information that do the definitions. We 
might have arrived at these simple expressions by arguing that they are sensible approximations in our context. These 
relations might be the result of our simplifying the information contained in A and B. With these last equations we 



can integrate equation (C4| over u> and fully get rid of any w-dependence, namely 

d t a(t,x,p)+v p -d x a(t,x,p) =/ 8(e p - e pl )r pp ,a(t,x,p') . (C7) 

J P ' 

This is a closed eqution for a without any explicit appearance of A and B or of other w dependent quantities. Without 



something like the expressions (C6| we would have not been able to arrive at a closed equation such as (C7) 
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